clear
close all
tic
%% note: DEFAULT UNIT: mm, ns, deg
%% pi has already defined by matlab

%% constant parameters of the flat detector
pitch = 0.4;
x_strip_width = 0.26; % the width of the x strips
y_strip_width = 0.35; % the width of the x strips
width = x_strip_width*1.5; % the 1.5 is the range of weight field
gaswidth = 5; % width of the gas gap

%% constant parameters of the gas
% % ar:co2 = 85:15
% sigma_trans = 0.191; % mm/sqrt(cm)
% sigma_longi = 0.172;
% velocity = 0.02336; % mm/ns
% lorentz_angle = 29.3; % deg

% % ar:ic4h10 = 90:10
% sigma_trans = 0.336; % mm/sqrt(cm)
% sigma_longi = 0.198;
% velocity = 0.0465; % mm/ns
% lorentz_angle = 29.3; % deg

% ar:CO2:CF4 = 45:15:40 % 560V/cm, according to huangquxuan
sigma_trans = 0.127; % mm/sqrt(cm)
sigma_longi = 0.123;
velocity = 0.0224; % mm/ns
lorentz_angle = 22.4; % deg
%% constant parameters of the simulation
step_x_pos = 0.025; %the step of the XPOS prob. and other distributions
x_min = -5; % lower edge of the detector in x axis
x_max = 5; % upper edge of the detector in x axis
N_x_pos = (x_max - x_min)/step_x_pos + 1;

step_z_pos = 0.1; % drift length divide
z_min = 0;%0;
z_max = gaswidth;
N_z_pos = (z_max - z_min)/step_z_pos + 1;

step_alpha = 5; % deg divide
alpha_min = 0;
alpha_max = 0;
N_alpha = (alpha_max - alpha_min)/step_alpha + 1;

step_drift_time = 1; % drift time divide
delta_time = step_drift_time/2;
drift_time_min = -20;
drift_time_max = round(gaswidth/cos(lorentz_angle/180*pi)/velocity) + 30 ;
N_drift_time = (drift_time_max - drift_time_min)/step_drift_time + 1;

strip_electron_min= 1;
strip_electron_max = 200;
strip_electron_array = [strip_electron_min; 2; 3; 5; 10; 20; 50; 100; 150; strip_electron_max]; 
N_strip_electron = size(strip_electron_array,1);

x0 = 0;
z0 = 4.5;
alpha = 10; % track incident angle %%note, cannot be 0!!!!
x_error_threhsold = 0.01; % the threshold in the normalized prob hist, to define the range as the x error
z_error_threhsold = 0.01; % the threshold in the normalized prob hist, to define the range as the x error
time_bin_width = 1; % the bin width in time calculation, especitally for the time advance correction.
prob_polya_over_threshold = 1.; % prob of polya distribution over amp threshold
%% calculation: the prob. distribution of x pos (influenced by transverse diffusion)
% key target distributions or functions:
x_error = zeros(N_z_pos, N_alpha);
x_error_mean = zeros(N_z_pos, 1);
z_error = zeros(N_z_pos, N_alpha);
time_advance_bias_temp = zeros(N_strip_electron, N_z_pos);
time_advance_bias = zeros(N_strip_electron, N_alpha);

temp_time_line = zeros(N_drift_time, 1);

% here we make a standard gaussian distribution, and get the cdf of it.
std_gaussian_min = -5;
std_gaussian_max = 5;
step_gaussian = 0.001;
standard_gaussian = makedist('Normal','mu',0,'sigma',1);
std_gaussian_array = std_gaussian_min : step_gaussian : std_gaussian_max;
cdf_std_gaussian= cdf(standard_gaussian,std_gaussian_array);
% plot(std_gaussian_array,cdf_std_gaussian)

x_pos_array = (x_min : step_x_pos : x_max)'; 
z_pos_array = (z_min : step_z_pos : z_max)';  %record the z pos array
alpha_array = (alpha_min : step_alpha : alpha_max)';  %record the alpha array
drift_time_array = (drift_time_min : step_drift_time : drift_time_max)';  %record the drift_time array


%alpha_array(1)=2;
prob_lower = zeros(N_x_pos,N_z_pos);
prob_upper = zeros(N_x_pos,N_z_pos);
prob_reaching = zeros(N_x_pos,N_z_pos); % prob. of a electron arrives the weight field range
prob_reaching_normalized =  zeros(N_x_pos,N_z_pos); % normalized prob. 

% first, for different z, the longitudinal diffusion caused time
% distribution is various. get the distribution in each z.
drift_time_z_pos = zeros(N_drift_time, N_z_pos);
for j = 1:1:N_z_pos
    z0 = z_pos_array(j);
    z0cor = z0/cos(lorentz_angle/180*pi);
    time0 = z0cor/velocity;
    for k = 1:1:N_drift_time
        %lower edge of drfit time 
        tempt = (drift_time_array(k) - time0 - delta_time)/(sigma_longi * sqrt(z0cor/10) / velocity);
        if ( tempt >= std_gaussian_min && tempt <= std_gaussian_max)
            drift_time_lower = cdf_std_gaussian(round((tempt - std_gaussian_min)/step_gaussian) + 1);
        else
            drift_time_lower = 0;
        end
        %upper edge of drift time 
        tempt = (drift_time_array(k) - time0 + delta_time)/(sigma_longi * sqrt(z0cor/10) / velocity);
        if ( tempt >= std_gaussian_min && tempt <= std_gaussian_max)
            drift_time_upper = cdf_std_gaussian(round((tempt - std_gaussian_min)/step_gaussian) + 1);
        else
            drift_time_lower = 0;
            drift_time_upper = 0;
        end
        drift_time_z_pos(k,j) = drift_time_upper - drift_time_lower;
        
    end
end

% calculate the x error and z error for a certain alpha, various z
for i = 1 : 1 : N_alpha
    % refresh alpha
    alpha = alpha_array(i);
    
    % calculate the ± time difference in a strip, with alpha
    if alpha  + lorentz_angle ~= 0
        strip_time_half_difference = abs(pitch / 2 * cos(lorentz_angle/180*pi) / tan((alpha + lorentz_angle)/180*pi) / velocity);
    else
        strip_time_half_difference = abs(z_max/2 / cos(lorentz_angle/180*pi) / velocity);
    end
    
    if (abs(alpha + lorentz_angle) >= 8)
        for j = 1 : 1 : N_z_pos
            %refresh z0
            z0 = z_pos_array(j);
            z0cor = z0/cos(lorentz_angle/180*pi);
            for k = 1 : 1 : N_x_pos  
                if z0cor + (x_pos_array(k) - x0) * cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi) >= 0 && z0cor + (x_pos_array(k) - x0) * cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi) <= gaswidth/cos(lorentz_angle/180*pi)
                    %calculate the prob.lower
                    tempx = (x_pos_array(k) - x0 - width/2)*cos(lorentz_angle/180*pi)/(sigma_trans * sqrt((z0cor + (x_pos_array(k) - x0) * cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi))/10));
                    if (tempx >= std_gaussian_min && tempx <= std_gaussian_max)
                        prob_lower(k,j) = cdf_std_gaussian(round((tempx - std_gaussian_min)/step_gaussian) + 1);
                    else
                        if (tempx < std_gaussian_min)
                            prob_lower(k,j) = 0;
                        else
                            prob_lower(k,j) = 1;
                        end
                    end
                    %calculate the prob.upper
                    tempx = (x_pos_array(k) - x0 + width/2)*cos(lorentz_angle/180*pi)/(sigma_trans * sqrt((z0cor + (x_pos_array(k) - x0) * cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi))/10));
                    if (tempx >= std_gaussian_min && tempx <= std_gaussian_max)
                        prob_upper(k,j) = cdf_std_gaussian(round((tempx - std_gaussian_min)/step_gaussian) + 1);
                    else
                        if (tempx < std_gaussian_min)
                            prob_upper(k,j) = 0;
                        else
                            prob_upper(k,j) = 1;
                        end
                    end

                else
                    prob_lower(k,j) = 0;
                    prob_upper(k,j) = 0;
                    prob_reaching(k,j) = 0;
                end 
            end

        end

        prob_reaching =  prob_upper - prob_lower;

        sum_prob_reaching = sum(prob_reaching);
        prob_reaching_normalized = prob_reaching ./ sum_prob_reaching;

        % find the lower and upper pos of the prob over threshold and define the
        % x error range and z error range
        for j = 1 : 1: N_z_pos
            % for x error
            x_error_lower = x_min - 1;
            x_error_upper = x_min - 1;
            x_prob_accumulation = 0;
            % lower edge
            for k = 1 : 1 : N_x_pos
                x_prob_accumulation = x_prob_accumulation + prob_reaching_normalized(k,j);
                if (x_prob_accumulation >= x_error_threhsold)
                    x_error_lower = x_pos_array(k);
                    break;
                end
            end
            %upper edge
            x_prob_accumulation = 0;
            for k = N_x_pos : -1 : 1
                x_prob_accumulation = x_prob_accumulation + prob_reaching_normalized(k,j);
                if (x_prob_accumulation >= x_error_threhsold)
                    x_error_upper = x_pos_array(k);
                    break;
                end
            end

            if (x_error_lower > x_min && x_error_upper > x_min)
                x_error(j,i) = max(abs(x_error_upper - x0), abs(x_error_lower - x0));% - abs(abs(x_error_upper) - abs(x_error_lower))/2;%
            end
            
            % for z error
            % the prob_reaching_normalized is also used for z error
            temp_time_line = zeros(N_drift_time, 1);
            z0 = z_pos_array(j);
            z0cor = z0/cos(lorentz_angle/180*pi);
            % focus on the prob. distribution within the range
            for k = find(x_pos_array == x_error_lower) : 1 : find(x_pos_array == x_error_upper)
                if x_pos_array(k) - x0 <= 0
                    tempz = z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi); %mm
                else
                    tempz = z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi) + (x_pos_array(k) - x0)*sin(lorentz_angle/180*pi); %mm
                end
                if round((tempz * cos(lorentz_angle/180*pi) - z_min)/step_z_pos) + 1 <= N_z_pos && round((tempz * cos(lorentz_angle/180*pi) - z_min)/step_z_pos) + 1 >= 1
                    temp_time_line = temp_time_line + prob_reaching_normalized(k,j) * drift_time_z_pos(:,round((tempz * cos(lorentz_angle/180*pi) - z_min)/step_z_pos) + 1); % prob*hist
                else
                    if round((tempz * cos(lorentz_angle/180*pi) - z_min)/step_z_pos) + 1 <= N_z_pos
                        temp_time_line = temp_time_line + prob_reaching_normalized(k,j) * drift_time_z_pos(:,N_z_pos); % prob*hist
                    else
                        temp_time_line = temp_time_line + prob_reaching_normalized(k,j) * drift_time_z_pos(:,1); % prob*hist
                    end
                end
            end
            
            % now, the time distribution of a certain z, certain alpha is
            % calculated.
            % for z error
            drift_time_lower = drift_time_min - 1;
            drift_time_upper = drift_time_min - 1;
            drift_time_prob_accumulation = 0;
            % lower edge
            for k = 1 : 1 : N_drift_time
                drift_time_prob_accumulation = drift_time_prob_accumulation + temp_time_line(k);
                if (drift_time_prob_accumulation >= z_error_threhsold)
                    drift_time_lower = drift_time_array(k);
                    break;
                end
            end
            %upper edge
            drift_time_prob_accumulation = 0;
            for k = N_drift_time : -1 : 1
                drift_time_prob_accumulation = drift_time_prob_accumulation + temp_time_line(k);
                if (drift_time_prob_accumulation >= z_error_threhsold)
                    drift_time_upper = drift_time_array(k);
                    break;
                end
            end

            if (drift_time_lower > drift_time_min && drift_time_upper > drift_time_min)
                drift_time_error = (drift_time_upper - drift_time_lower)/2; % ns
                % max(abs(drift_time_upper - z0/velocity), abs(drift_time_lower - z0/velocity));%
                % then convert to z error (mm)
                z_error(j,i) = drift_time_error * velocity;
                %z_error(j,i) = max(abs(x_error_upper), abs(x_error_lower));% - abs(abs(x_error_upper) - abs(x_error_lower))/2;%
            end
            
            % for time advance correction
            % the time advance distribution is set: x0 , t=0 (in a certain strip, with mid z0)
            temp_time_line_2 = (-(ceil(strip_time_half_difference) + 101) : 1 : (ceil(strip_time_half_difference) + 101))';
            prob_once_in_t = zeros(2 * ceil(strip_time_half_difference) + 203, 1);
            prob_another_slower_than_t = zeros(2 * ceil(strip_time_half_difference) + 203, 1);
            
            z0 = z_pos_array(j);
            z0cor = z0/cos(lorentz_angle/180*pi);
            % focus on the prob. distribution within the range
            for k = find(x_pos_array == x_error_lower) : 1 : find(x_pos_array == x_error_upper)
                % get the temp z for a x pos
                if x_pos_array(k) - x0 <= 0
                    tempz = z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi); %mm
                else
                    tempz = z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi) + (x_pos_array(k) - x0)*sin(lorentz_angle/180*pi); %mm
                end
                
                % now, the alpha, z0, and zk is determined. so in the residual
                % distribution prob_once_in_t and
                % prob_another_slower_than_t,
                % the t=0 (no bias time) is set as the REAL T in Z0. 
                real_time = z0cor/velocity;
                % think about the longitudinal diffusion (±3sigma) in zk:
                time_bias = (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi)/velocity;
                longi_difference = 3*sigma_longi*sqrt(tempz/10)/velocity;
                
                for n = 1 : 1 : 2 * ceil(strip_time_half_difference) + 203
                    % calculate the prob of once hit in time t, and the
                    % prob of another hit slower than t
                    if (temp_time_line_2(n) - time_bin_width/2 - time_bias)/(longi_difference/3) >= std_gaussian_min && (temp_time_line_2(n) + time_bin_width/2 - time_bias)/(longi_difference/3) <= std_gaussian_max
                        prob_once_in_t(n) = prob_once_in_t(n) + prob_reaching_normalized(k,j) * (cdf_std_gaussian(round(((temp_time_line_2(n) - time_bias + time_bin_width/2)/(longi_difference/3) - std_gaussian_min)/step_gaussian) + 1) - cdf_std_gaussian(round(((temp_time_line_2(n) - time_bias - time_bin_width/2)/(longi_difference/3) - std_gaussian_min)/step_gaussian) + 1));                                        
                        prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + prob_reaching_normalized(k,j) * (1 - cdf_std_gaussian(round(((temp_time_line_2(n) - time_bias)/(longi_difference/3) - std_gaussian_min)/step_gaussian) + 1));
%                         prob_once_in_t(n) = prob_once_in_t(n) + prob_reaching_normalized(k,j) * (cdf_std_gaussian(round(((temp_time_line_2(n) - m + time_bin_width/2)/(sigma_longi*sqrt(tempz/10)/velocity) - std_gaussian_min)/step_gaussian) + 1) - cdf_std_gaussian(round(((temp_time_line_2(n) - m - time_bin_width/2)/(sigma_longi*sqrt(tempz/10)/velocity) - std_gaussian_min)/step_gaussian) + 1));              
%                         prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + prob_reaching_normalized(k,j) * (1 - cdf_std_gaussian(round(((temp_time_line_2(n) - m)/(sigma_longi*sqrt(tempz/10)/velocity) - std_gaussian_min)/step_gaussian) + 1));
                    else
                        prob_once_in_t(n) = prob_once_in_t(n) + 0;
                        if (temp_time_line_2(n) - time_bias - time_bin_width/2)/(longi_difference/3) < std_gaussian_min
                            prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + prob_reaching_normalized(k,j) * 1;
                        else
                            prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + 0;
                        end
                    end    
                end
                 
            end
            if max(prob_once_in_t) > 0
                prob_once_in_t = prob_once_in_t / sum(prob_once_in_t);
            end
            prob_another_slower_than_t = prob_another_slower_than_t / max(prob_another_slower_than_t);
            
            % calculate the time advance bias, with a certain alpha, z0,
            % and various Nelectron
            prob_temp_time_advance = zeros(2 * ceil(strip_time_half_difference) + 203, 1);
            for k = 1 : 1 : N_strip_electron
                prob_temp_time_advance = strip_electron_array(k) * prob_polya_over_threshold .* prob_once_in_t .* (prob_another_slower_than_t .* prob_polya_over_threshold + (1 - prob_polya_over_threshold)) .^(strip_electron_array(k)-1);
                temp_weight_sum = sum(prob_temp_time_advance .* temp_time_line_2);
                temp_weight = sum(prob_temp_time_advance);
                time_advance_bias_temp(k,j) = round(temp_weight_sum/temp_weight);
                % this is wrong! when z is small, and N is small, the
                % prob_temp_time_advance is NOT GAUSSAIN !!!
                % temp_pos_array = find (prob_temp_time_advance == max(prob_temp_time_advance));
                % time_advance_bias_temp(k,j) = temp_time_line_2(temp_pos_array((1+size(find (prob_temp_time_advance == max(prob_temp_time_advance)),1))/2),1);
            end
            
        end
        
        
    else
        % if the alpha + lorentz_angle == 0, so the x error equals the pitch/2
        for j = 1 : 1: N_z_pos
            x_error(j,i) = pitch/2;
            z_error(j,i) = 1/2;
        end
        
        for j = 1 : 1 : N_z_pos
            % for time advance correction
            % the time advance distribution is set: x0 , t=0 (in a certain strip, with mid z0)
            temp_time_line_2 = (-(ceil(strip_time_half_difference) + 101) : 1 : (ceil(strip_time_half_difference) + 101))';
            prob_once_in_t = zeros(2 * ceil(strip_time_half_difference) + 203, 1);
            prob_another_slower_than_t = zeros(2 * ceil(strip_time_half_difference) + 203, 1);
            
            z0 = z_pos_array(j);
            z0cor = z0/cos(lorentz_angle/180*pi);
            % focus on the prob. distribution within the range
            for k = find(x_pos_array == 0) - (pitch/2/step_x_pos) : 1 : find(x_pos_array == 0) + (pitch/2/step_x_pos)
                % get the temp z for a x pos
%                 if x_pos_array(k) - x0 >= 0
%                     tempz = (z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi)); %mm
%                 else
%                     tempz = 0;
%                 end
                
                if x_pos_array(k) - x0 <= 0
                    tempz = z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi); %mm
                else
                    tempz = z0cor + (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi) + (x_pos_array(k) - x0)*sin(lorentz_angle/180*pi); %mm
                end
                if tempz < z_min/cos(lorentz_angle/180*pi)
                    tempz = 0;
                end
                if tempz > z_max/cos(lorentz_angle/180*pi)
                    tempz = z_max/cos(lorentz_angle/180*pi);
                end
                
                % now, the alpha, z0, and zk is determined. so in the residual
                % distribution prob_once_in_t and
                % prob_another_slower_than_t,
                % the t=0 (no bias time) is set as the REAL T in Z0. 
                real_time = z0cor/velocity;
                % think about the longitudinal diffusion (±3sigma) in zk:
                time_bias = (x_pos_array(k) - x0)*cos(lorentz_angle/180*pi)/tan((alpha + lorentz_angle)/180*pi)/velocity;
                longi_difference = 3*sigma_longi*sqrt(tempz/10)/velocity;
                
                for n = 1 : 1 : 2 * ceil(strip_time_half_difference) + 23
                    % calculate the prob of once hit in time t, and the
                    % prob of another hit slower than t
                    if (temp_time_line_2(n) - time_bin_width/2 - time_bias)/(longi_difference/3) >= std_gaussian_min && (temp_time_line_2(n) + time_bin_width/2 - time_bias)/(longi_difference/3) <= std_gaussian_max
                        prob_once_in_t(n) = prob_once_in_t(n) + (cdf_std_gaussian(round(((temp_time_line_2(n) - time_bias + time_bin_width/2)/(longi_difference/3) - std_gaussian_min)/step_gaussian) + 1) - cdf_std_gaussian(round(((temp_time_line_2(n) - time_bias - time_bin_width/2)/(longi_difference/3) - std_gaussian_min)/step_gaussian) + 1));                                        
                        prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + (1 - cdf_std_gaussian(round(((temp_time_line_2(n) - time_bias)/(longi_difference/3) - std_gaussian_min)/step_gaussian) + 1));
%                         prob_once_in_t(n) = prob_once_in_t(n) + prob_reaching_normalized(k,j) * (cdf_std_gaussian(round(((temp_time_line_2(n) - m + time_bin_width/2)/(sigma_longi*sqrt(tempz/10)/velocity) - std_gaussian_min)/step_gaussian) + 1) - cdf_std_gaussian(round(((temp_time_line_2(n) - m - time_bin_width/2)/(sigma_longi*sqrt(tempz/10)/velocity) - std_gaussian_min)/step_gaussian) + 1));              
%                         prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + prob_reaching_normalized(k,j) * (1 - cdf_std_gaussian(round(((temp_time_line_2(n) - m)/(sigma_longi*sqrt(tempz/10)/velocity) - std_gaussian_min)/step_gaussian) + 1));
                    else
                        prob_once_in_t(n) = prob_once_in_t(n) + 0;
                        if (temp_time_line_2(n) - time_bias - time_bin_width/2)/(longi_difference/3) < std_gaussian_min
                            prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + 1;
                        else
                            prob_another_slower_than_t(n) = prob_another_slower_than_t(n) + 0;
                        end
                    end                     
                end
                 
            end
            if max(prob_once_in_t) > 0
                prob_once_in_t = prob_once_in_t / sum(prob_once_in_t);
            end
            
            prob_another_slower_than_t = prob_another_slower_than_t / max(prob_another_slower_than_t);
            
            % calculate the time advance bias, with a certain alpha, z0,
            % and various Nelectron
            prob_temp_time_advance = zeros(2 * ceil(strip_time_half_difference) + 203, 1);
            for k = 1 : 1 : N_strip_electron
                prob_temp_time_advance = prob_polya_over_threshold .* prob_once_in_t .* (prob_another_slower_than_t .* prob_polya_over_threshold + (1 - prob_polya_over_threshold)) .^(strip_electron_array(k)-1);
                temp_weight_sum = sum(prob_temp_time_advance .* temp_time_line_2);
                temp_weight = sum(prob_temp_time_advance);
                time_advance_bias_temp(k,j) = round(temp_weight_sum/temp_weight);
                % temp_pos_array = find (prob_temp_time_advance == max(prob_temp_time_advance));
                % time_advance_bias_temp(k,j) = temp_time_line_2(temp_pos_array((1+size(find (prob_temp_time_advance == max(prob_temp_time_advance)),1))/2),1);
            end
        end
    end
    
    % now, the time_advance_bias_temp is obtained, and it should be fitted
    % in 2-D, and get some parameters..?
    % or, using a 3D hist to store the data, and make a 3D fitting....
    % much more complex
%     figure (i)
%     plot(strip_electron_array, time_advance_bias_temp);
end

x_error_mean = mean(x_error')';



figure (1)
plot(z_pos_array, x_error)
figure (2)
plot(z_pos_array, x_error_mean)
figure (3)
plot(z_pos_array, z_error)
figure (4)
plot(strip_electron_array, time_advance_bias_temp);
% 
% figure (5)
% plot(temp_time_line_2, prob_once_in_t, 'r');
% 
% figure (6)
% plot(temp_time_line_2, prob_another_slower_than_t, 'b');
% figure (7)
% plot(temp_time_line_2, prob_temp_time_advance, 'r');

% now, fit the x error by self-defined function
x_error_func = @(a,b,c,x) a*x.^2 + b*x + c;
x_error_cfg = fittype(x_error_func); % define the fit type
x_error_init_guess = [-0.004, 0.05, 0.1]; %define the init of the 3 paras
x_error_cf= fit(z_pos_array, x_error, x_error_cfg, 'StartPoint', x_error_init_guess);
x_error_fit_result = coeffvalues(x_error_cf)


% temp_zposarray = zeros(N_z_pos, N_alpha);
% temp_alphaarray = zeros(N_z_pos, N_alpha);
% for i = 1 : 1 : N_alpha
%     temp_zposarray(:,i) = z_pos_array;
% end
% for i = 1 : 1 : N_z_pos
%     temp_alphaarray(i,:) = alpha_array';
% end

% now, the z error fitting is not complete yet...
% z_error_cfg = fittype(@(a, b, c, d, e, x, y) (a*x.^2 + b*x + c).*(d*y.^e), 'independent', {'x', 'y'}, 'dependent', 'z'); % define the fit type
% 
% % g = fittype( @(a,b,c,d,x,y) a*x.^2+b*x+c*exp(-(y-d).^2), ...
% %             'independent', {'x', 'y'}, ...
% %             'dependent', 'z' ); % for fitting surfaces
% z_error_init_guess = [-0.1, 0.6, 0.4, 10, -1]; %define the init of the 3 paras        
% z_error_cf= fit([temp_zposarray temp_alphaarray], z_error, z_error_cfg, 'StartPoint', z_error_init_guess);
% z_error_fit_result = coeffvalues(z_error_cf);

z_error_func = @(a,b,c,x) a*x.^2 + b*x + c;
z_error_cfg = fittype(z_error_func); % define the fit type
z_error_init_guess = [-0.004, 0.05, 0.1]; %define the init of the 3 paras
z_error_cf= fit(z_pos_array, z_error, z_error_cfg, 'StartPoint', z_error_init_guess);
z_error_fit_result = coeffvalues(z_error_cf)

time_advance_test(:,1) = -time_advance_bias_temp(:,2);
time_advance_test(:,2) = -time_advance_bias_temp(:,6);
time_advance_test(:,3) = -time_advance_bias_temp(:,11);
time_advance_test(:,4) = -time_advance_bias_temp(:,16);
time_advance_test(:,5) = -time_advance_bias_temp(:,21);
time_advance_test(:,6) = -time_advance_bias_temp(:,26);
time_advance_test(:,7) = -time_advance_bias_temp(:,31);
time_advance_test(:,8) = -time_advance_bias_temp(:,36);
time_advance_test(:,9) = -time_advance_bias_temp(:,41);
time_advance_test(:,10) = -time_advance_bias_temp(:,46);
time_advance_test(:,11) = -time_advance_bias_temp(:,51);